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Abstract 

We present an extended version of THERMINATOR, a Monte Carlo event generator 
dedicated to studies of the statistical production of particles in relativistic heavy- 
ion collisions. The increased functionality of the code contains the following features: 
The input of any shape of the freeze-out hypersurface and the expansion velocity 
field, including the 3+1 dimensional profiles, in particular those generated externally 
with various hydrodynamic codes. The hypersufraces may have variable thermal 
parameters, which allows for studies departing significantly from the mid-rapidity 
region, where the baryon chemical potential becomes large. We include a library of 
standard sets of hypersurfaces and velocity profiles describing the RHIC Au+Au 
data at ^Jsnn = 200 GeV for various centralities, as well as those anticipated for the 
LHC Pb+Pb collisions at y/sNN = 5.5 TeV. A separate code, FEMTO-THERMINATOR, 
is provided to carry out the analysis of femtoscopic correlations which are an im- 
portant source of information concerning the size and expansion of the system. We 
also include several useful scripts that carry out auxiliary tasks, such as obtain- 
ing an estimate of the number of elastic collisions after the freeze-out, counting of 
particles flowing back into the fireball and violating causality (typically very few), 
or visualizing various results: the particle pr-spectra, the elliptic flow coefficients, 
and the HBT correlation radii. We also investigate the problem of the back-flow of 
particles into the hydrodynamic region, as well as estimate the elastic rescattering 
in terms of trajectory crossings. The package is written in C++ and uses the CERN 
ROOT environment. 



Key words: relativistic heavy-ion collisions, statistical hadronization, Monte Carlo 
event generator, LHC, RHIC, SPS, FAIR, NICA 
PACS: 25.75.-q, 25.75.Dw, 25.75.Ld 



Preprint submitted to Elsevier 



Jan. 30, 2011 



Program summary 



Title of the program: THERMINATOR 2 

Catalog identifier: 

Program web page: 

|ht tp : / / therminator 2 . if j . edu . pi/ 

Licensing provisions: none 

Computer: any computer with a C++ compiler and the CERN ROOT environ- 
ment, ver. 5.26 or later, tested with Intel Core2 Duo CPU E8400 @ 3 GHz, 
4 GB RAM 

Operating system under which the program has been tested: 
Linux Ubuntu 10.10 x64 (gcc 4.4.5) ROOT ver. 5.26 

Programming language used: C++ with the CERN ROOT libraries, BASH shell 
No. of lines in distributed program, including test data: 
13000 code, 300 configuration, 400000 input files (hypersurface library) 
No. of bytes in distributed program, including test data: 2.7 MB 
Distribution format: tar.gz 

Nature of physical problem: particle production via statistical hadronization 
in relativistic heavy-ion collisions 

Method of solution: Monte-Carlo simulation, analyzed with ROOT 
Restrictions concerning the complexity of the problem: none 
Memory required to execute with typical data: 

- 30 MB therm2_events, 

- 150 MB therm2_f emto 

Typical running time: default configuration at 3 GHz 

- primordial multiplicities 70 min (calculated only once per case), 

— 8 min/500 events, 

— 10 min - draw all figures, 

- 25 min/one kr bin in the HBT analysis with 5000 events 
Typical data file size: default configuration 

— 45 MB/500 events 

- 35 MB/correlation file (one fcy bin) 

- 45 kB/fit file (projections and fits) 



* Supported in part by the Polish Ministry of Science and Higher Education, grant 
No. N N202 263438 and Foundation for Polish Science. 

Email addresses: Mikolaj . ChojnackiSif j .edu.pl (Mikolaj Chojnacki), 
kisiel@if.pw.edu.pl (Adam Kisiel), Wojciech.Florkowski@ifj.edu.pl 
(Wojciech Florkowski), Wojciech.Broniowski@if j .edu.pl (Wojciech 
Broniowski) . 
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1 Introduction 



We introduce an updated and largely extended version of THERMINATOR PQ, the 
THERMal heavy IoN generATOR, created to carry out the statistical hadroniza- 
tion in relativistic heavy-ion collisions. Numerous successful analyses have 
been performed with the help of our code over the last few years [2H63]. con- 
tributing to better understanding of the involved evolution of the hot and 
dense system created in ultra-relativistic heavy-ion collision. 

THERMINATOR 2 is a Monte Carlo generator written in C++ and using the stan- 
dard CERN ROOT [M] environment. That way, apart from model applications, 
the code can be easily adapted for purposes directly linked to experimental 
data analysis, detector modeling, or estimates for the heavy-ion experiments 
at RHIC, LHC, SPS, FAIR, or NICA. 

Originally, THERMINATOR was designed to perform the hadronic freeze-out on 
simple boost-invariant hypersurfaces, such as the Cracow single-freeze-out 
model [HS1ES] and the Blast- Wave model [dTHd^] . Meanwhile, the code has 
evolved into a versatile tool, where the freeze-out profile and the expansion 
velocity field of any shape can be implemented, allowing application to all ap- 
proaches based on statistical hadronization on a specified hypersurface. The 
present functionality of the code includes the following new features: 

(1) An implementation of any shape of the freeze-out hypersurface and the 
expansion velocity field is possible now, including the 2+1 and 3+1 di- 
mensional profiles, in particular those generated externally with various 
codes for perfect hydrodynamics [22|.[T0l48T] . 

(2) The hypersufraces may have space-time dependent thermal parameters, 
which allows for studies of non-boost-invariant systems. In particular, one 
may depart significantly from the mid-rapidity region to the fragmenta- 
tion regions where the baryon chemical potential becomes large [8||4"5]. 

(3) The package includes a library of "standard" sets of hypersurfaces and 
velocity profiles, which describe the Au+Au data [22] at the highest RHIC 
energy ^snn = 200 GeV for various centralities. We also provide the 
hypersurfaces and velocity profiles anticipated [23] for the LHC Pb+Pb 
collisions at v /J/v r /v = 5.5 TeV, which prior to studies with the real data 
can be used for the modeling of detectors. 

(4) A separate code, the FEMTO-THERMINATOR, is provided to carry out the 
analysis of the femtoscopic correlations. These correlations are an impor- 
tant source of complementary information concerning the size and expan- 
sion of the system, which should be reproduced in a realistic description 
of the data [5j[22j|26]. 

We also provide several scripts that carry out auxiliary tasks: 
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(1) Counting the number of particles flowing back into the fireball, thus vio- 
lating causality (for realistic freeze-out profiles this number is negligible, 
below 1%). 

(2) An estimate of the number of elastic collisions (trajectory crossings) after 
the freeze-out. For realistic cases this number is small (~ 1.5 for pions 
scattering off other pions), supporting the single- freeze-out approxima- 
tion. 

(3) Visualization of the results: the particle p^-spectra, the elliptic flow co- 
efficient t>2, or the HBT correlation radii. 

We have also made an effort to make the package more user-friendly, providing 
numerous examples of usage in typical situations. 

Other codes implementing the thermal description of hadronization have also 
been developed recently, to mention SHARE |S21IH3] and THERMUS [HI] , comput- 
ing the particle abundances, DRAGON [2B] , a Monte Carlo generator of particle 
production from a fragmented fireball, or HYDJET++ [53 | l85" | l86]. a heavy- ion 
event generator with hydrodynamics and jets. A code similar to THERMINATOR, 
Fast Hadron Freeze-out Generator, has been presented in Refs. [8T|[88] . 



2 Statistical hadronization 

For the detailed description of the physics behind THERMINATOR 2 the reader 
is referred to the original paper [TJ. Here we only wish to stress that the 
statistical approach to heavy-ion collisions has proved to be very useful for a 
very broad range of observables. With few physical input parameters, such as 
the temperature, chemical potentials, size, and the velocity of the collective 
flow, the models describe the observed particle abundances [65 , 68 , 69 , 1891 - 1107] . 
the transverse-momentum spectra [66], the balance functions |108] ll09j. the 
elliptic flow coefficient |110]llllj . and femtoscopic observables [5l[22j|26]. A 
uniform description of the spectra, the elliptic flow coefficient t>2, and the 
HBT data achieved in [22J explained for the first time the so-called "RHIC 
HBT puzzle". 

The key element of the success behind the statistical approach is the inclusion 
of the full list [112] of hadronic resonances, whose number grows rapidly ac- 
cording to the Hagedorn hypothesis [113H118] . At the rather high temperature 
of the freeze-out, ~ 140—165 MeV, the resonances contribute very significantly 
to the observables. THERMINATOR 2 uses the particle data tables in the univer- 
sal form originally included in the SHARE [82J package, and incorporates the 
four-**** and three-*** resonances. 

All particles are created at the freeze-out hypersurface according to the Cooper- 
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Frye formula |119j . The two- and three-body decays of resonances, proceeding 
in cascades, produce the stable particles observed in the detectors. THERMINATOR 
offers the full space-time information on positions and velocities of the pro- 
duced particles. As a Monte Carlo event generator written in C++ in the CERN 
ROOT [M] environment, THERMINATOR can be straightforwardly interfaced to 
the standard software routinely used in the data analysis of relativistic heavy- 
ion collisions at the LHC, RHIC, SPS, FAIR or NICA. The inclusion of exper- 
imental acceptance and kinematic cuts poses no difficulty in a Monte Carlo 
generator. 



2. 1 Cooper-Frye formalism 



Our starting point for the analysis of different models is the famous Cooper- 
Frye formula [119] . It expresses the number of hadrons being produced on the 
freeze-out hypersurface S M by the following integral 

N=(2s + l)f J d^(x)jff(x,p), (1) 

where 



f(p ■ u) 



exp 



p^ - (B/jb + J 3 /i /3 + Sfj. s + Cfx 



c 



T 



± 1 



(2) 



is the phase-space distribution function of the particles (the stable ones and 
resonances), and d£ M is the three-dimensional element of the freeze-out hy- 
persurface. The latter may be calculated with the help of the formula known 
from the differential geometry [120] 



d£ M = e m py-^-Qp--g—dad(3dy, 



(3) 



where e Ma/ 3 7 is the Levi-Civita tensor (£0123 = +1) and the variables a, (3, 
and 7 are used to parametrize the three-dimensional freeze-out hypersurface 
submerged in the four-dimensional Minkowski space. The quantity dH^ may 
be interpreted as the four-vector normal to the hypersurface with the norm 
equal to the "volume of the hypersurface element". Explicitly, 



oE 



dx dx dx 

da 80 

dy dy dy 

da d/3 df 

dz dz dz 

da 8/3 df 



dadfidy, 



(4) 
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Fig. 1. The system of coordinates used to parametrize the boost-invariant freeze-out 
hypersurfaces obtained from the 2+1 hydrodynamic codes. The dashed line shows 
the freeze-out points with the fixed values of the azimuthal angle <j>. 

while the remaining components are obtained via cyclic permutations of t, z, 
y, and z. In the following subsections we shall explicitly provide dS M in the 
form appropriate for various considered cases. The formulas will be directly 
used to generate the hadrons with the Monte-Carlo method. 

The quantities appearing in the expressions below have the following physical 
interpretation: 

• thermodynamic quantities on the freeze-out hypersurface: 

• T - temperature, 

• Hb - baryon chemical potential, 

• /i/ 3 - isospin chemical potential, 

• fis - strange chemical potential, 

• fie - charmed chemical potential, 

• space-time coordinates: 

■ r = \Jt 2 — z 2 - longitudinal proper time, 

• p = \Jx 2 + y 2 - distance in the transverse plane, 

■ <p - azimuthal angle, 

■ Y s — ~ ln((i + z)/(t — z)) - space-time rapidity, 

• coordinates and parameters used to specify the freeze-out hypersurface: 

■ n - the initial proper time for hydrodynamics, 

■ Tf - the final proper time, 

■ ( - angle in p — r plane, see Fig. [T] 

• 6 - angle between the Yg-axis and the direction determined by the space 
time point (rj, 0,0,0) and the point on the freeze-out hypersurface, see 
Fig.El 

■ d = d((,(j),9) - distance from the space-time point (r,, 0,0,0) to the point 
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on the freeze-out hypersurface, see Fig. [2j 

• flow characteristics of the fluid element: 

• u x , u y - transverse components of the four- velocity, 

• vt = \Jv 2 + Vy - magnitude of the transverse velocity, 

• (f)f - azimuthal angle of the transverse velocity, 

• Yf - rapidity of the fluid element, 

• properties of a particle generated on the freeze-out hypersurface: 

• s - spin, 

• m - mass, 

• Pt - transverse momentum, 

• rriT = J m 2 + p\ - transverse mass, 

■ <f) p - azimuthal angle of the transverse momentum, 

■ Y - rapidity, 



In the case of freeze-out hypersurfaces obtained in the framework of the 
perfect-fluid hydrodynamics the distribution function in (CQ) is the Fermi-Dirac 
or Bose-Einstein distribution and the number of hadrons produced is given by 
the formula 



(0 , ; dY p T dp T d<p p f dZ^p 



{2s + 1) J J 



(2vr)3 J exp(- bvb+Is^+Shs+Cm, ) exp(^) ± 1 ' 

(5) 



T-Ti 




Fig. 2. The system of coordinates used to parametrize the freeze-out hypersurfaces 
obtained from the 3+1 hydrodynamic codes, see Eqs. ([6]) and ([?]). 
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2.2 Freeze-out in 3+1 -dimensional perfect- fluid hydrodynamics 



First, we consider the most general case where both the freeze-out hypersurface 
and the hydrodynamic flow at freeze-out are obtained from the 3+1 perfect- 
fluid hydrodynamic code, with no symmetries assumed. In this situation the 
freeze-out hypersurface is parametrized in the following way |35j, 



+ ( i A<r a a\ ■ a ■ r\ i, d{C„<j>,0) cos 6 
t—[Ti + d[Q, 0, 9) smt/sin Q cosh — , 

x = d((, (j), 9) sin 9 cos ( cos 4>, 
y = d((, <f>, 9) sin 9 cos ( sin cf>, 

z=( Ti + d((, 0, 9) sin 9 sin C) sinh d ^'^^ COs6 , ( 6 ) 



where d is a function of the parameters 0, (, and 9, which is obtained numer- 
ically from the hydrodynamic code, see Fig. [2j The parameter A is a scale 
used to change the dimensionless space-time rapidity Y s into a dimensional 
quantity (for example one may choose A = Pmax/Y™ 8 '*), namely 



Y s = Y s {(, 0, 9) = , 

T = r(C, 0, 9) = n + d(C, 0, 9) sin 9 sin (. (7) 



In addition, we express the particle four-momentum and the fluid four-velocity 
in terms of particle and fluid rapidities, respectively, which leads to the ex- 
pression 



P^ = Jl + u 2 x + u 2 [m T cosh(Y/ — Y)— Pt(u x cos <p p + u y sin (f> p )] . (8) 



The calculation of the integration measure in the Cooper-Frye gives 
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dE^ = d 2 sinfl 



(r» + dsin£) 
A 



dd ( ( dcosO" 

— cos C — rn T cos C cosh Y — 

d( \ V A . 



+Pt sin ( cos ((p — 4> p 



'a ■ a dd a 
d sin 6 — — cos o 



( , , ( dcos9\ > / 

x rriT sin £ cosh Y — + px cos Q cos( 



A 



( r n ^ . „\ . „ . , /„ dcos#' 

+ d cos & + — sin & sin sinn r : — 

V 89 V A 



(9) 



Using (J7J) in ([2]) one obtains a more compact form 



d sin 9 
dE^jT = — d t 



dd 



cos ( (—rriT cos ( cosh (Y — Y s 



dd 



+p T sin C cos (0 — 4> p )) + cos £ sin 9 yi sin — — cos 
x (niT sin C cosh (F — Y s ) + pr cos £ cos (0 — 



H — [ d cos 9 + ^ sin # | sin # sinh (y — Y s 

r 89 Vs. 



(10) 



Freeze-out in 2+1 -dimensional boost-invariant hydrodynamics 



For boost-invariant systems, considered in approximate treatment of the heavy- 
ion collisions, the fluid rapidity is equal to space-time rapidity, Yf = Y s , and 
the function d depends only on and (. In this case we find [18J 



\-v% 



(rriT cosh(F — Y s ) — v T p T cos(0 — (f> p )) 



(11) 



and 
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T 



freeze-out curve 

T = Tf + A p 



Pmax 



Fig. 3. The freeze-out curves considered in the modified Blast- Wave model [5]. 



dd 

— cos ( (—m T cos ( cosh (Y — Y s ) + p T sin ( cos (0 — P )) 



+rf cos £ (ray sin £ cosh (Y — Y" s ) + cos ( cos (0 — <j> p )) 
+ 7T7 Pt sm(0 - 0, 



(12) 



2.4 Blast-wave models 



In this Section we discuss one of the most popular parameterizations of the 
freeze-out hypersurface — the Blast- Wave model. In its standard form, the 
model is boost-invariant and cylindrically symmetric. Moreover, the assump- 
tion is made that the freeze-out happens at a constant value of the proper 
time 

r (C) = T f — const. (13) 
In order to get a broader applicability, we generalize this condition to the 
formula 

r(p) = Tf + Ap, (14) 

where A is a constant that describes the slope of the freeze-out curve in the 
Minkowski space, see Fig. [3J With A > (A < 0) we may consider the freeze- 
out scenarios where the outer parts of the system freeze-out later (earlier). Of 
course, with A = we reproduce the standard Blast-Wave parametrization. 
In this case we use 



--{rriT cosh(F s — Y) — vtPt cos(0 — (f) p )) 



(15) 
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Pmax 



P 

Fig. 4. The freeze-out curve assumed in the Cracow model. 

and 

= (t/ + Ap)p (m T cosh(y s — F) — Ap T cos(0 - <j> p )) . (16) 

The freeze-out conditions defined by Eq. ffl4|) with different values of A were 
studied in [5]. 

2.5 Cracow model 

Another boost-invariant and cylindrically symmetric freeze-out model [66] as- 
sumes that the hypersurface is defined by the condition, 

afx^ = t 2 -x 2 -y 2 - z 2 = t 2 - p 2 = t 2 = const., (17) 

see Fig.Hl and the fluid four- velocity is proportional to the space-time position 

u" = 7(l,v) = - = -(l (18) 

Tf Tf \ t) 

Note that the physical interpretation of the parameter r/ is now slightly dif- 
ferent from that used in the previous Section. 

The calculation of the volume element of the freeze-out hypersurface shows 
that it is proportional to u^, as in the model of Siemens and Rasmussen [121] . 
so we may write 

= — (yrj + p 2 m T cosh(F s -Y) - ppr cos(0 - P )) (19) 
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and 



dE^ = r f ppX- (20) 
2.6 Resonance decays 

The two- and three-body resonance decays are implemented in exactly the 
same manner as in the original version of the code pp. The format of the 
tables from SHARE [H2] is presented in Appendix [A] 



3 Structure of THERMINATOR 2 

THERMINATOR is written in C++ and uses the CERN ROOT environment - An 
object-oriented data analysis framework [M] . 

3. 1 Modules and their functions 

The code consists of several modules performing subsequent tasks. First, the 
information on particles and their decays is read by Parcer . cxx and stored 
by ParticleDB . cxx. Then EventGenerator . cxx conducts the generation of 
events, which consists of several tasks. Integrator . cxx computes the aver- 
age multiplicity of primordial particles corresponding to a given model and 
parameters. Next, Event. cxx generates the event, first producing the pri- 
mordial particles with the multiplicity given by the Poisson distribution from 
Integrator . cxx, and then calling ParticleDacayer . cxx to carry out the 
decays of resonances. Finally, the event is written to the output file. The cor- 
responding block structure is shown in Fig. [5j The organization of the input 
and output files is shown in Fig. |6j 

3.2 File structure 

After unpacking the THERMINATOR 2 package in the directory therminator2 
one obtains the directory structure presented in FigJTl The main directory of 
THERMINATOR 2 contains the following files: 

Doxygen - Doxygen |122] configuration, 

events.ini - THERMINATOR 2 configuration used by therm2_events, 
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Configurator 

contains therm2_events settings 



Parser 

reads configuration file 



ParticleDB 

ParticleType database containing DecayTable 
with two and three particle DeeayChannels 



Parser 

reads SHARE files 



! share/particles. data i 
| share/decays.data ! 



EventGenerator 

generates desired number of events 
and saves them to files 



Integrator 

calculates average particle type multiplicities 
provides primordial particles 



| even'ts/*/fmuliticity_<hash>~txt 



Model (abstract) 

creates a particle on a model 
dependet freezeout hypersurface 



' fomodel/'.ini 1 
i fomodel/lhyquid'/'.xml • 



Event 

contains a list of Particles 



ParticleDecayer 

decay of resonances 



i events/'/event.txt 1 
! events/VeventOOO.root ' 
! events/'/event001 .root i 

Fig. 5. Visual concept of the therm2_events program. Solid line boxes represent 
classes or a groups of classes with their names in bold font. Dashed rectangular 
boxes represent input or output files. The arrows symbolize the communication 
between classes. 



femto.ini - THERMINATOR 2 configuration used by therm2_f emto, 
hbtfit.ini - THERMINATOR 2 configuration used by therm2_hbtf it, 
Makefile - make configuration, 
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SHARE data tables 

share/particles.data 
share/decays.data 



configuration file 
events.ini 



ROOT TTree type event files 

events/*/eventOOO.root 
events/*/event001 .root 



configuration file 
femto.ini 








ROOT histograms 
events/*/femto*.root 







configuration file 
hbtfit.ini 





freeze-out model parameters 

fomodeir.ini 

fomodel/lhyquid3d/*.xml 

fomodel/lhyquid2dbi/*.xml 



auxiliary file 
events/7fmultiplicity_*.txt 



ASCII event file (optional) 
events/7event.txt 



physical results, figures 

events/7fig_*.eps 


i 






FEMTO-THERMINATOR 



Fig. 6. Organization of input and output files. Rectangular boxes represent files, 
ovals correspond to codes or scripts. The arrows indicate reading or writing. The 
part enclosed in the dashed line corresponds to the FEMTO-THERMINATOR package. 

manual . html - reference manual in HTML format (created with make doc), 
redirection to doc/html/index.html, 

manual.pdf - reference manual in PDF format (created with make doc), 
symbolic link to doc/latex/refman.pdf, 
runtorque.sh - bash shell script for TORQUE [123J, 
runall . sh - bash shell script for the full calculation, 
runevents . sh - bash shell script for generating events, 
runf igure . sh - bash shell script for generating figures, 
runhbt . sh - bash shell script for the femtoscopic analysis, 
therm2_events - program generating events (created with make), 
therm2_f emto - program for two-particle correlations (created with make), 
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-therminator2 / 
-addons/ 
-build/ 

-doxygen/ 

-include/ 

-obj/ 

Lsrc/ 
-doc/ 
-events/ 
-f omodel/ 

l-lhyquid2dbi/ 

Llhyquid3d/ 
-macro/ 
-share/ 

Fig. 7. Directory structure of the THERMINATOR 2 package. 
therm2_hbtf it - program for HBT analysis (created with make). 

The other directories are: 

therminator2/build/doxygen/ 

- contains Doxygen input files used to generate the reference manual, 
therminator2/build/include/ 

- contains header files (*.h) with the declarations of classes and structures, 
as well as the general declaration file THGlobal.h, 

therminator2/build/obj / 

- contains the compiled source files (* . o) and the binaries therm2_events and 
therm2_f emto (created with make), 

therminator2/build/src/ 

- contains the source code files for the classes (*.cxx) and the programs 
therm2_events . cxx and therm2_f emto . cxx, 

therminator2/doc/ 

- contains the Doxygen reference manuals in the HTML and LaTeX formats 
(created with make doc), 

therminator2/events/ 

- contains the generated event files, figures, correlation function files, etc. (the 
contents of the events/ directory and the subdirectory structure depends on 
the selected model and its settings). 

The default subdirectories for each model are as follows: 

therminator2/event/krakow/ - the Krakow Single Freeze-Out model, 
therminator2/event/blastwave/ - the Blast- Wave model, 
therminator2/event/bwa/ - the Blast- Wave A-class model, 
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therminator2/event/lhyquid2dbi-*/ - the Lhyquid 2+1D Boost-Invariant 
hydro, 

therminator2/event/lhyquid3d/rn - the Lhyquid 3+1D hydro. 
The files stored in the events/*/ subdirectory: 
event*. root - event file, 

f multiplicity.* .txt - primordial particle multiplicity, 

f ig_* . eps - figure in EPS format, 

f ig_* .png - figure in PNG format, 

f ig_* .xml - figure data in XML format, 

f ig_* . C - ROOT macro to generate a figure, 

femto*.root - two particle correlation function, 

hbtf it* . root - fit to the two particle correlation function, 

events.ini - event generator configuration file, 

femto.ini - two-particle correlation configuration file, 

hbtf it . ini - HBT fit configuration file, 

model_config.ini - freeze-out model configuration file (the file name de- 
pends on the chosen model, e.g. Ihyquid2dbi . ini). 

File naming conventions are as follows: The event*. root files are always 
stored in a sequence beginning with eventOOO . root, and subsequently with 
the file counter incremented by one. The f multiplicity.* .txt file has an 
eight-character hex-base number representing the CRC-32 identification num- 
ber calculated from the model parameters, e.g., f multiplicity_10D57D490 . txt 
for the default case of RHIC with the centrality 20-30%. The f ig_* . * files have 
the name inherited from the corresponding ROOT macros (with the word "fig- 
ure" shortened to "fig"). The femto*.root and fmfit*.root filenames de- 
pend on the particle pair type, average pair momentum, and the inclusion of 
particles form resonances, e.g., f emtopipi0200a. root for the pion-pion pairs 
with average kx = 200 MeV, and with resonance decays idicated with "a" . 

therminator2/f omodel/ 

- contains the freeze-out model configuration files (* . ini) and freeze-out hy- 
persurfaces data in the XML format (*.xml). The five freeze-out model con- 
figuration files are supplied: 

blastwave . ini - for Blast- Wave model, 
bwa. ini - for modified Blast- Wave A-Class model, 
krakow . ini - for Krakow Single Freeze-out model, 
lhyquid2dbi . ini - for 2+1 boost-invariant hydrodynamics, 
lhyquid3d.ini - for 3+1 hydrodynamics. 



THERMINATOR 2 distribution does not include hypersurfaces for this case. They 
can be downloaded from the project homepage http : //therminator2 . if j . edu.pl 
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The freeze-out hypersurface input files for RHIC at ^Jsnn = 200 GeV and 
various centralities have the following names: 

lhyquid2dbi/RHICAuAu200c0005Ti500ti025Tf 145. xml (for c = 0-5 %), 
lhyquid2dbi/RHICAuAu200c0510Ti491ti025Tf 145. xml (5-10 %), 
lhyquld2dbi/RHICAuAu200cl020T1476ti025Tf 145. xml (10-20 %), 
lhyquid2dbi/RHICAuAu200c2030Ti455ti025Tf 145. xml (20-30 %), 
lhyquid2dbi/RHICAuAu200c3040Ti429ti025Tf 145. xml (30-40 %), 
lhyquid2dbi/RHICAuAu200c4050Ti398ti025Tf 145. xml (40-50 %), 
lhyquid2dbi/RHICAuAu200c5060Ti354ti025Tf 145. xml (50-60 %), 
lhyquid2dbi/RHICAuAu200c6070Ti279ti025Tf 145. xml (60-70 %). 

The name encodes the colliding nuclei, the collision energy (in MeV), the 
centrality class, the initial central temperature Ti (in MeV), the starting time 
used in hydrodynamics ti (in fm/c), and the freeze-out temperature Tf (in 
MeV). In an analogous way, the freeze-out hypersurface input files for the 
LHC collisions at ^snn = 5.5 TeV and various centralities are: 

lhyquid2dbi/LHCPbPb5500c0005Ti500til00Tf 145. xml (for c = 0-5 %), 

lhyquid2dbi/LHCPbPb5500c0510Ti500til00Tf 145. xml (5-10 %), 

lhyquid2dbi/LHCPbPb5500cl020Ti500til00Tf 145. xml (10-20 %), 

lhyquid2dbi/LHCPbPb5500c2030Ti500til00Tf 145. xml (20-30 %), 

lhyquid2dbi/LHCPbPb5500c3040T1500til00Tfl45.xml (30-40 %), 

lhyquid2dbi/LHCPbPb5500c4050T1500til00Tfl45.xml (40-50 %), 

lhyquid2dbi/LHCPbPb5500c5060T1500til00Tfl45.xml (50-60 %), 

lhyquid2dbi/LHCPbPb5500c6070T1500til00Tfl45.xml (60-70 %). 

therminator2/macro/ 

- contains the ROOT macro files for generating figures in the event analysis: 

f igure_disteta. C - distribution in pseudorapidity, 

f igure_distmt . C - distribution in transverse mass, 

f igure_distpt . C - distribution in transverse momentum, 

f igure_distpt_exotic . C - distribution in transverse momentum of less 

abundant particles (e.g., those with strangeness), 

f igure_distpt_pion.C - pions' anatomy in transverse momentum (contri- 
butions from individual resonance decays to pion spectra), 
f igure_distpt_proton . C - protons' anatomy in transverse momentum (con- 
tributions from individual resonance decays to proton spectra), 
f igure_di strap . C - distribution in rapidity, 

f igure_di strap _prim. C - primordial-particle distribution in rapidity, 

f igureJisret . C - primordial-particle distribution in (p,t),t), 

f igure_hsxyt . C - primordial-particle distribution in (x,y,r), 

f igure_xpx.C - position-momentum correlation, x — p x , 

f igure_v2pt . C - differential elliptic flow at midrapidity, 

f igure_v4pt . C - differential V4 coefficient at midrapidity, 

f igure_v2rap . C - integrated elliptic flow vs. rapidity. 

Other ROOT macros used by the figure*. C files are: 

events2chain.C - reads event files, 

hist2xml.C - converts histogram data into the XML format, 
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model21egend.C - creates figure legends displaying the model parameters. 

therminator2/share 

- contains the SHARE [22] database files: 

decays. data - particle decay channels, 
particles .data - particle properties. 

therminator2/addons 

- contains an auxiliary program (estimate of the elastic rescattering) : 
coll .mk - makefile, 

collcount . cxx - program estimating elastic collisions per pion. 
3.3 *.ini files 

The main configuration file events . ini (Table [1]) is used by therm2_events 
program to select the model, number of generated events, type of output, file 
and directory names, etc. Depending on the value of FreezeOutModel, the 
program reads additional information from an appropriate model settings file 
* . ini. The file is located in the FreezeOutDir directory. 

The model configuration file (Table [2]) contains the parameters of a partic- 
ular model, e.g., the configuration file for the default FreezeOutModel = 
Lhyquid2DBI is f omodel/lhyquid2dbi . ini. 

The variable IntegrateSamples, used to obtain the integrand ranges and the 
particle average multiplicities, is set to a large value (5000000 by default). 
This leads to accurate results, but increases the running time of the first stage 
of therm2_events. The user may decrease this value, however, care is needed 
in the inspection of results. A very large value of NumberOf Events, above 
20000, is needed to achieve stable HBT results. A low value (about 1000) is 
sufficient to obtain the pr spectra. Of course, the running time is proportional 
to NumberOf Events. 



4 Storage of events 

The particles are stored in the computer memory as collections corresponding 
to one event. Once the whole event is generated (all primordial particles are 
produced and all unstable particles have been sequentially decayed), the event 
is written down to the output file and erased from the memory. There are two 
types of output: plain ASCII text and ROOT file. 



18 



4-1 Plain text output file 



The plain text output has the advantage of having the complete information 
and being easily readable by a human-being or a computer program. It is, 
however, not efficient in terms of disk-space usage and data-access time. The 
structure of the plain text file is the following: 

• The first two lines form the file header containing a description of the 
EVENT_ENTRY, see further in text. 

• Each event in the output file consists of an event header containing the 
EVENT_ID (a unique number generated with the CRC-32 algorithm from 
the system time and event sequence number) and number of particles in the 
event. 

• The lines that follow the event header are the particle entries. Each particle 
is represented by a single line. 

An example of a few lines of the test output: 

• THERMINATOR 2 text output 

#<EVENT_ENTRY>eid fathereid pid fatherpid rootpid decayed mass e px py pz t x y z</EVENT_ENTRY> 

#<EVENT_ID>0x65679F6</EVENT_ID> 

#<N0_0F_PARTICLES>9684</NQ_0F_PARTICLES> 

-1 9001 9001 9001 2.350000e+00 5.423351e+00 5.105272e-01 8.400604e-01 -4 . 787891e+00 

3.012708e+01 3.797786e+00 1.937673e+00 -2 . 789553e+01 

1 -1 5218 5218 5218 1 2.250000e+00 5.035699e+00 1 . 177745e+00 -3 . 007674e-01 -4 . 337997e+00 

2.272015e+01 2.773765e+00 -1 . 740429e+00 -1 . 985397e+01 

The particle event entry columns contain the following information: 

• eid, the particle sequence number in the event, 

• fathereid, the sequence number of the parent, 

• pid, Particle Data Group (PDG) identification number, 

• fatherpid, PDG identification number of the parent, 

• rootpid, PDG identification number of the original primordial particle 
starting the decay sequence, 

• the decay flag (1 if particle has decayed, if not), 

• mass and components of the four momentum: e, px, py, and pz (in GeV), 

• space-time coordinates of the creation point: t, x, y and z (in fm/c). 

The information on eid, fathereid, together with pid and fathereid can 
be used to trace back the resonance cascade from the last stable particle up 
to its primordial particle born on the freeze-out hypersurface. Particles with 
fathereid = -1 are primordial particles with pid = fathereid = rootpid. 



19 



4-2 ROOT file 



The information stored in the ROOT file is slightly different from the one 
stored in the text file and it is grouped in three separate TTree type objects. 

The tree called particles consists of one branch named particle. The branch 
is filled with particle entries (leaves) - one leaf per particle. This tree stores 
information on particles from 500 events (by default). The particle leaf consist 
of a class called ParticleCoor that stores the following particle properties: 



data type 


name 


description 


Float.t 


mass 


mass 


Float.t 


t, x, y, z 


space-time coordinates (in fm/c) 


Float.t 


e, px, py, pz 


energy and momentum coordinates (in GeV) 


Int.t 


decayed 


decay flag 


Int.t 


pid 


PDG identification number 


Int.t 


fathcrpid 


parent PDG number 


Int.t 


rootpid 


root (primordial) particle PDG number 


Int.t 


cid 


sequence number in the event 


Int.t 


fathereid 


parent sequence number in the event 


Ulnt.t 


cvcntid 


CRC-32 unique id of the event 



The data types used here are the ROOT data types. 

The events tree called events consists of one branch named event. The branch 
is filled with event entries (leaves) - one leaf per event, by default 500 entries 

per file. The event leaf consist of a structure called StructEvent that stores 
the information on the number of particles in each event and the event id 
number: 

data type name description 
Ulnt.t eventID unique ID of the event generated with CRC-32 algorithm 

Ulnt.t entries number of particle entries in this event 

Ulnt.t entriesprev number of previous entries, by default. 

The entriesprev variable is filled with and is not used at this stage. It is 
used in the classes that read a sequence of event files for analysis. 

The parameters tree called parameters consists of eight branches that store 
the following information : 



20 



branch name 


data type 


description 


IntegrateSamplc 


Ulnt.t 


number of Monte-Carlo integration samples 


Randomize 


Ulnt.t 


random seed switch (1-random seed, 0-constant seed) 


TimcStamp 


Char.t 


date and time of this run 


ModellD 


Ulnt.t 


id number of the chosen freeze-out model 


ModelName 


Char.t 


name of the freeze-out model 


ModelHash 


Char.t 


CRC-32 hash made of the freeze-out model parameters 


ModclDcscription 


Char.t 


model parameter description 


ModclParameters 


ModeLt 


structure containing model parameters. 



The structure ModeLt is the default freeze-out model structure that does 
not contain any parameter information. Each freeze-out model has a separate 

structure defined and saved to the file. The detailed information can be found 
in the therminator2/build/include/StructModel .h file. As an example, 

the structure for the default freeze-out model is presented as follows: 



data type 


name 


description 


Float.t 


RapPRange 


rapidity range 


Float.t 


RapSRangc 


spatial rapidity range 


Float.t 


Taul 


initial proper time 


Float.t 


TempF 


freeze-out temperature 


Float.t 


MuB 


baryochemical potential 


Float.t 


Mul 


isospin potential 


Float.t 


MuS 


strangeness potential 


Float.t 


MuC 


charm potential 


Float.t 


CollidingEnergy 


CMS beam energy 


Float.t 


CentralityMin 


lower limit of ccntrality 


Float.t 


CentralityMax 


higher limit of centrality 


Float.t 


ImpactParameter 


impact parameter 


Float.t 


Tbmpl 


initial central temperature 


Char.t 


DeviceName [30] 


name of the device, e.g. RHIC, LHC 


Char.t 


CollidingSystem[30] 


colliding nuclei, e.g, Au, Pb. 



5 Installation and sample runs 

5.1 Downloading the latest version 

THERMINATOR 2 is distributed in a form of a .tar.gz archive containing the C++ 
sources, scripts, configuration files, and the data files. The latest package can 
be downloaded from the website 

http : //therminator2 . if j . edu .pi/ . 

Prior to compiling and running THERMINATOR 2, the user has to install the 
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following programs: 



C++ compiler, 

ROOT ver. 5.26 or later. 



The latest version of the ROOT package is available at http : / / root . cern . ch/ 
The simplest way to proceed is to follow the instructions below. In this way 
one first creates a directory called therminator2. Then, the package is saved 
and extracted, 

mkdir therminator2 
cd therminator2/ 

wget therminator2 . if j . edu.pl/therminator2-latest . tar . gz 
tar xzf therminator2-latest . tar . gz 

In place of "latest" the actual number appears. The information about the 
version of the package is contained in the version file. 



5.2 Compilation 

One compiles THERMINATOR 2 with the command 
make 

After a successful compilation, the binary files are present in the main direc- 
tory. 



5.3 Doxy gen documentation 



The source files of THERMINATOR 2 are documented with Doxygen to generate 
the web-page and the PDF files with detailed information on the classes, 
variables, and files. The user must have the Doxygen package installed on 
the system (latest version can be obtained from www.doxygen.org/). The 
documentation is created with the command 

make doc 

In the main distribution directory, two files are created: manual.html, which 
redirects the web browser to the main documentation page, and manual.pdf, 
which is a symbolic link to the PDF reference manual. 
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5.4 A test run 



A test run is executed with with the script 
. /runall . sh 

This run displays all basic functionality of the package, including the femto- 
scopic analysis, producing the plots with the physical results, etc. It corre- 
sponds to Au+Au collisions at the top RHIC energy and centrality 20-30%. 
The 2D Boost-Invariant hydrodynamics is used as the input. By default, 5000 
events are used. The test run takes several hours to complete (4 h 20 min on 
a 3 GHz processor). At the end the events and figures are contained in the 
events/lhyquid2dbi-RHICAuAu200c2030Ti455ti025Tf 145/ directory. 

5. 5 Running 

Now we describe in a greater detail the elements (codes and figure macros) 
used in the script of the previous subsection. 

The main code of THERMINAT0R 2 is therm2_events, which generates the 
events. It has been set by default (see the description of the * . ini files) 
to generate 500 events for the case of the RHIC top energy and centrality 
20-30%. The 2D Boost-Invariant hydrodynamics is used. In order to perform 
the default run, issue the command 

. /therm2_events 

This creates the file eventnnn . root in the subdirectory 
therminator2/events/lhyquid2dbi-RHICAuAu200c2030Ti455ti025Tf 145 
together with a text file fmultiplicity_10D57D490.txt that holds the pri- 
mordial particle multiplicities. The number identifier nnn is 000 if no event 
files are present in the subdirectory, whereas it is one unit larger from the file 
with highest identifier if the files are already present. 

The script 

. /rune vents . sh 

generates as many events as requested by the Number Of Events parameter in 
the events.ini file. The user may add more events by repeating the . /run- 
events, sh command. 

To generate via a ROOT macro the figure showing the pt spectra, issue the 
command 
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root -x ' . /macro/f igure_distpt . C(" . /events/ 

lhyquid2dbi-RHICAuAu200c2030Ti455ti025Tf 145/" ,n) ' 

where n is the number of event files used, starting from eventOOO . root. The 
single and double quotation marks are necessary to pass the parameters to the 
macro. The macro creates two files in the events subdirectory: f ig_distpt . eps 
and f ig_distpt .xml. The first one is the EPS graphics file, the second one 
contains the data from the histogram in a form of the XML file. 

The creation of all figures from ROOT macros can be done automatically by 
issuing the command 

. /runf igure . sh . /events/lhyquid2dbi-RHICAuAu200c2030Ti455ti025Tf 145/ 

All executables, i.e. binaries and shell scripts, produce help information when 
called with parameter -h or — help, for instance, 

. /therm2_events -h 
./runf igure — help 



5.6 Auxiliary files 

After the program therm2_events is started, the file event_* . tmp is created in 
the main directory of the THERMINATOR 2 distribution. This file has two lines 
with the information about the subdirectory with the event files generated by 
the current run and the total number of event files located in that directory, 
including files from the previous runs. The name of the temporary file ends 
with the shell parent process ID (PPID), i.e. the PID number of the shell script 
that has started the therm2_events process. The auxiliary event_*.tmp file 
is used by scripts in order to be able to run multiple THERMINATOR 2 processes 
at once, e.g., on a computer cluster. 



6 Scripts to generate physical results 

We provide several ROOT scripts that may be used to generate physical results. 
These scripts read the chain of event* . root files and produce figures and ROOT 
figure macros. To work properly, the scripts must read from the sequence of 
event files generated by THERMINATOR 2, starting from eventOOO . root. This 
means that the user should not accidentally delete event files from the middle 
of a sequence. 
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The transverse- momentum spectra are obtained by running the figure_sp.C 
script in ROOT. To use n event files, the user executes the following command 
(from the therminator2 directory): 

root -x ' . /macro/f igure_distpt . C(" . /events/ 

lhyquid2dbi-RHICAuAu200c2030Ti455ti025Tf 145/" ,n) ' 

where n is the requested number of the event files to be processed. 

For the elliptic flow coefficient, t>2, the user should run 

root -x ' . /macro/f igure_v2pt . C(\" . /events/ 

lhyquid2dbi-RHICAuAu200c2030Ti455ti025Tf 145/" ,n) ' 

This computes the elliptic flow coefficient, t>2, as a function of the transverse 
momentum, px, in the reaction plane set by the hydrodynamics. The script 
for the fourth-order Fourier flow coefficient, t>4, called f igure_v4pt . C, is run 
analogously. 



7 HBT analysis 



This section describes the new functionality of THERMINATOR 2, allowing to 
compute the HBT correlation radii with the two-particle method. 
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Fig. 8. Transverse-momentum spectra of ir + , K + , and protons for Au+Au collisions 
at y/sNN = 200 GeV and the centrality class 20-30%, fig_distpt.eps (protons 
from A decays are excluded). 
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Fig. 9. The elliptic flow coefficient, v 2 , for tt + , K + , and protons at the RHIC top 
energy ^snn = 200 GeV and the centrality class 20-30%, fig_v2pt.eps. 

7.1 Method 



The method of carrying out the femtoscopic analysis is the same as described 
in detail in Ref. 0. For integrity of this paper, we present below the main 
points. 

Consider the two-particle distribution expressed by the two-particle emission 
function, 



W 2 (pi,P2) = E Pl E P2 — J S(x 1 ,x 2 ,pi,p 2 )d*x 1 d 4 x 2 . 

The correlation function is then defined as 

W 2 (pi,p 2 ) 



C(pi,p 2 ) 



Wi(pi)Wi(p 2 



where 



Wl (p) = E p ^- = Jd 4 xS(x }P ) } 
with S(x,p) being the source emission function. 
We define the momentum difference of the particles, 



(21) 



(22) 



(23) 



(24) 
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the sum of their momenta, 

P = (P , P) = (E P1 + E^n + p 2 ) , (25) 
and the average momentum of the pair, 

k = (fc , k) = l - (E P1 + E^fr+fa) . (26) 
For the particles with equal masses we use the notation 

<finv = 2k*. (27) 

The space and time separations of the members of the pair are r = r \ — r 2 
and At = t% — ti- When calculated in the Pair Rest Frame (PRF), they are 
denoted as r* and At*. Both k* and f* appear as the arguments of the wave 
function, since PRF is the natural reference frame for the representation of 
the wave-function structure. 

In general, the HBT analysis may be performed in any reference frame. One 
determines the correlation function as a function of the relative momentum 
components in the selected frame. Then the inverse widths of the correlation 
functions yields the size parameters of the system in this frame. Here we use 
the standard Bertsch-Pratt decomposition [12411125] of the mean and relative 
three-momenta into the three components. The long axis coincides with the 
beam axis, the out axis is determined by the direction of the average transverse 
momentum of the pair, denoted later by kx, and the side direction is perpen- 
dicular to the other two axes. Following the RHIC experiments, we choose 
to perform the analysis in the longitudinal co-moving system (LCMS), which 
is defined as a system where fci ong = 0. In this work the notation is used in 
which the values in PRF are denoted by an asterisk, while the values without 
asterisk are defined in LCMS. 

By definition, in a Monte-Carlo method the numerical equivalent of the inte- 
grals (I21|23p is the summation over particles or pairs of particles generated 
by the Monte-Carlo procedure. The numerical calculation of the correlation 
functions is done in bins, which may be expressed with the help of the function 



Sa(p) 



1 if \Px\ < f , \Py\ < f , \Pz\ < T 



(2* 



otherwise. 



Then the correlation function may be expressed as 
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C(q, k) 



E EM? -Pi +Pj)S A (k - Hpi+p^ik*, 



i j^i 



E E S A (q - Pi + Pj)S A (k -\{pi + Pj)) 



(29) 



Various analyses of the HBT correlations use different approximations for 
the full pionic wave function, In the non- interacting system, or in the 
interacting but non-relativistic case, the motion of the center of mass can be 
separated and one deals with the relative motion only. The simplest relative 
wave function ignores all dynamical interactions and has the form 

*=_L( e ^ + e -**^) j (30) 

where symmetrization over the two identical particles has been performed. 
Therefore |#| 2 = 1 + cos (2jfeV). 

Correlation functions calculated according to (129]) represent the ideal Bose- 
Einstein correlation functions. They are also very useful in model studies, 
because they can be calculated analytically for simple Gaussian emission func- 
tions. 

In addition, in the applied method mixing of particles is introduced when 
evaluating the correlation functions. The purpose of generating mixed events 
is to break all correlations other than those originating from femtoscopy, for 
instance the correlations due to resonances. The details are explained care- 
fully in Sect. V of [48]. For the case of identical charged particles this mixing 
technique is innocuous, as there are no resonances producing two identical 
same-charge particles. In addition, the mixing method increases the statistics, 
which is highly desirable in our simulations. The calculation without mixing 
is not implemented in the present code. 

The inclusion of Coulomb effects is described in [5]. 



7.2 Calculation of correlation functions extraction of the HBT radii 



The correlation functions are obtained through a numerical implementation 
of Eqs. f l2"$|) and (13"U|) (no Coulomb effects). As described in the previous sec- 
tions, particles generated by THERMINATOR 2 are grouped into events, as in 
experiment. In each event every charged pion is combined with every other 
pion of the same charge. For each pion pair, |\1/| 2 is calculated and added to 
the numerator of Eq. (1291) in a bin corresponding to the pair's g out , g S ide and 
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gi on g. At the same time, 1 is added to the denominator of Eq. (!29l in the 
corresponding bin. The resulting ratio yields the correlation function. 

By making a proper selection of single pions and pairs of pions one may study 
the correlation functions as functions of various variables. For instance, taking 
into account the pairs of particles within a certain total momentum range only, 
one obtains femtoscopic information for a given Jct- 

Usually, the single-particle emission function is postulated to be a static 3- 
dimensional ellipsoid with a Gaussian density profile 



Sm = iVexp (- Jf - - _ . (31) 

This source function is static, i.e., it does not depend on particle momentum. 
In this case the integral leads to the well known formula: 



C (k±, g ut, Qside, giong) = 1 + A exp -R 2 out (k ± )q i 



iout 



-^side(^J-)9side -^fong ( ^-L ) 9l 



ong 



(32) 



The quantities R on t, -Rside and -Ri on g, known as the "HBT radii", are the widths 
of the Gaussian approximation. It is important to emphasize that formula (13 2 j) 
is commonly used to fit the experimental data and to represent the results of 
the model calculations, although the experimental or model emission functions 
are frequently far from gaussians. In the simple fitting examples given in this 
paper we do not address this problem. More detailed discussion of such effects 
can be found in [5J. 



7.3 therm2_femto and therm2-fithbt programs 



The program that implements the two-particle femtoscopic formalism de- 
scribed above, see Eq. ( l29|) . is called therm2_f emto. The syntax of the com- 
mand is: 

./therm2_hbt <KTBIN> <EVENT_DIR> <EVENT_FILES> [FEMT0_INI] [PPID] 

where the parameter KTBIN=0 ,1,2,3 selects the transverse-momentum bin of 
the pair, EVENT_DIR is the directory where the event*. root files are stored, 
and EVENT_FILES is the number of the files to be taken. Parameter FEMTCLINI 
(optional) is the name configuration file. This parameter is by default is set to 
f emto . ini. The PPID parameter (optional) is the system's process ID number 
used by the shell scripts. By default it is equal to 0. 
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The kx bins correspond to the following ranges of the pair momentum: 

- (0.15, 0.25) GeV/c, 

1 - (0.25, 0.35) GeV/c, 

2 - (0.35, 0.45) GeV/c, 

3 - (0.45, 0.6) GeV/c. 

The output of the therm2_f emto code consists of the ROOT files f emto* .root, 
with histograms of the numerator and denominator, 

which can be used to calculate the correlation function. In addition we provide 
the projections of the pair emission function S on the three components of 
the pair relative separation r calculated in LCMS. They are saved as 2D 
histograms, with the magnitude of relative momentum qi nv at the other axis. 

The therm2_hbtf it program is an implementation of the fitting procedure, 
which uses the simplest fitting functional form ( 13"2"j) . The syntax of the com- 
mands is: 

therm2_hbtfit <C0RR_FILE> [HBTFIT.INI] [PPID] 

The first parameter C0RR_FILE is the file containing the numerator and de- 
nominator histograms. The second parameter (optional) is set by default to 
. /f ithbt . ini. It is the name of the configuration file, containing the fitting 
parameters and switches. The third parameter (optional) is the process ID 
number, used by shell scripts. 

The naming convention for the parameters is as follows: The label [*_f it] 
indicates a parameter to be fitted. The name of the parameter, e.g., Rout, is 
assigned its initial value from which the fit starts. The *FitType flag defines 
the type of fitting. If it is set to fixed, then this value is fixed to the starting 
value in the fitting procedure. If limit is given, the value can vary in the 
range given by the *Min and *Max parameters. Finally, if it is equal to free 
than the value can be fitted without constraints. 

The five parameters, Norm (set to 1), Lambda, Rout, Rside, and Rlong, cor- 
respond to the parameters of Eq. ( 1321) : the overall normalization, the A pa- 
rameter, and the three HBT radii (Rout, Rside, and Ri on g)- The parameter 
MaxFitRange contains the range of q values in which the fitting is performed. 

The output of the therm2_hbtf it program is a text file hbtradii .txt. It con- 
tains the values of the fitted parameters and statistical errors. Subsequent runs 
of the program append information to this file. The macro macro/ f igure_hbtradii . C 
reads that information and creates a plot of the HBT radii. The therm2_hbtf it 
program also produces the output ROOT file hbtf it* .root which should be in- 
spected in order to judge the quality of the fit. It contains the one- dimensional 
slices of the calculated correlation function, as well as the corresponding slices 
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of the fitted function. 

The histograms are named in the following way: 
<Type>_<Direction>_<Width> 

where Type can be CFproj for the correlation function slices or FITptoj for 
the fitted function slices. Direction is one of out, side or long, while Width 
can have a value of 1 - corresponding to a slice which has a width of 1 bin 
in the other two directions, 2 - a slice with the width of 4 bins and 3 - a 
slice with the width of 10 bins. The slices of the correlation function and 
the fitted function for the same direction and the same width should be 
plotted on top of each other in order to confirm that the fitting procedure 
correctly captured the width of the correlation. This is done by the macro 
macro/f igure_f itf emto . C. 

We warn the user that obtaining statistically stable results for the HBT radii 
requires sometimes a very large number of events, such as 20000 or even 50000. 
This requires an appropriate change of the parameter NumberOf Events in the 
events.ini file, leading to long execution times. 

The running sequence of the femtoscopic analysis has been automatized in the 
shell script . / runhbt . sh. The user should simply execute 

. / runhbt . sh 

which in addition generates the physical plots. An example is presented in 
Fig. [inland Fig. HU 



8 Yielding auxiliary results 

In this Section we present new analyses, possible with THERMINAT0R 2. In fact, 
they relate to fundamental physical aspects of the freeze-out model assumed 
in our approach. 

8. 1 Backflow 

The first effect concerns the so-called back-flow problem [126H128] . It origi- 
nates from the fact that for typical freeze-out hypersurfaces generated with 
hydrodynamics, particles generated with the Cooper-Frye formula may flow 
back into the hydrodynamics domain, thus violating causality. A typical rem- 
edy is the removal (used also in THERMINATOR 2) of the back-flowing particles 
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Fig. 10. Projections of the correlation function (symbols) calculated with the default 
settings for pions with kx in the (0.15, 0.25 GeV) bin. Au+Au, ^snn = 200 GeV, 
centrality class 20-30%, 50000 events (fig_cfprojpipiOa.eps). The lines show the 
fit projected the same way. 




0.2 0.3 0.4 0.5 0.2 0.3 0.4 0.5 

k T (GeV/c) k T (GeV/c) 



Fig. 11. The HBT radii and A as a function of the average pair kx calculated with 
the default settings, Au+Au, ^/sWn = 200 GeV, centrality class 20-30%, 50000 
mixed events (fig_hbtradii.eps). 

via the constraint p^dY,^ > 0, where d£ M is normal to the hyper surf ace. To 
see how big is the back-flow effect, we compute the back- flowing particles re- 
questing that p M <i£ M < 0, and compare the result to the case without this 
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Fig. 12. Same as Fig. [8] for the back-flowing particles, 5000 events, 
constraint, i.e. including all particles. 

To prepare the binary file of therm2_events generating the back-flowing parti- 
cles only, the user should run (this feature is implemented for the Lhyquid2dbi 
and Lhyquid3d models only) 

make all BACK_FL0W=1 

and then execute 
. /therm2_events 

The user should take care of the fact that the directory specified with EventSubDir 
does not include the fmultiplicity_* .txt file, such that the multiplicities 
are recomputed for the present case of the back-flowing particles. 

The result is shown in Fig. [121 Comparing with Fig. El we note that the back- 
flow effect is very small (less than 0.5% for the pions and even less for more 
massive particles) in our realistic calculation. This smallness is caused by large 
transverse velocity of the fluid, which inhibits the thermal motion of particles 
in the opposite direction |22j . 

8.2 Elastic rescattering 

The other considered effect concerns the elastic rescattering after freeze-out. 
In our single-freeze-out approximation it is assumed that the chemical and 
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thermal freeze-outs occur simultaneously. As a matter of fact, this is not the 
case, as the inelastic processes, responsible for the chemical equilibrium, have 
generally a smaller cross section than the elastic processes. In order to estimate 
how good this approximation is, one may consider a simple estimate of elastic 
rescattering based on trajectory crossings. 

Consider two particles, born at the freeze-out hypersurface at coordinates 
(ti,af() and (t 2 ,x%). Their trajectories are 



x 1 (t)=xl + (t-t 1 )v u x 2 (t) = x% + (t - t 2 )v 2 , (33) 
and the closest approach occurs at 



_ {Vl ~ V2) ■ (gg ~ g) f .V 

[Vi -V 2 ) ■ («! -V 2 ) 

Since this time cannot be earlier than t\ or t 2 , we define to = m ax(t m , t\, t 2 ). 
The distance squared at the closest approach is d 2 = (xi(t ) — x 2 (t )) 2 . The 
collision occurs if, geometrically, ird 2 < cr e \(s), where a e \(s) is the elastic cross 
section for the Mandelstam variable s. 

We have analyzed in detail the case of a charged pion scattering elastically 
off the other pions. The elastic cross section is computed with the help of the 
realistic parametrization of the pion-pion phase shifts from |129j . The result is 
that for most-central collisions the pion scatters elastically off the other pion 
on the average 1.8 times, while for the mid-central and peripheral collisions 
this value drops below 1. This shows that for the pions the effects of elastic 
rescattering after freeze-out are small and our approximation is appropriate. 

The user may carry out this analysis by running 

cd addons 

cp <dest_dir>/eventOOO.root . 
make -f coll.mk 
. /collcount 

where <destdir> denotes the directory with the previously generated event 
files for the model and centrality of interest. 

The code collcount . cxx may also serve as a simple example of an "off-line" 
analysis with the THERMINATOR 2 events. Its structures (e.g., concerning the 
reading of events) may be directly imitated in similar applications. 
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9 Conclusion 



We hope that THERMINATOR 2 will continue to be a useful and versatile tool 
for the heavy-ion community. In particular, the possibility of implementing 
the general freeze-out conditions (read from external files) allows the user to 
investigate in detail the predictions of the hydrodynamic approach followed 
with statistical hadronization (implemented in the single-freeze-out scenario). 
This approach was very successful in reproducing, in a uniform way, the RHIC 
data [22J and the LHC data at ^s NN = 2.76 TeV |13UH132| . It was also very 
helpful in a variety of analyses [2H63] . Therefore, it has proved to be realistic 
and may serve as a practical tool in other physical studies or in detector 
modeling. 

A new part of the package, FEMTO-THERMINATOR, carries out the femtoscopic 
analysis using the two-particle method. This is of practical importance, as 
complete analyses of the heavy-ion data typically consist of the studies of 
spectra, flow, and femtoscopy - all are now implemented in THERMINATOR 2. 

We note that a small modification of the approach presented in this paper, 
accounting for viscosity effects which draw the distribution functions away 
from the thermal equilibrium, may be included in the straightforward way. 
That way the code can be used as an afterburner in the increasingly popular 
hydrodynamic studies with viscosity effects [IT?lll34l4T4"8] . Fugacities, reflecting 
incomplete chemical equilibrium [133] of the system, can now be also incor- 
porated. This allows for model studies of this phenomenon. 

We have made a considerable effort to make the structure of the code clear and 
straightforward to adapt for other applications. The program is documented 
in detail, also via the Doxygen utility. Numerous ROOT macros, allowing to 
carry out various tasks such as the calculation and presentation of physical 
results, are supplied with the package. 

We foresee the following further extensions and applications of the program: 

(1) Statistical hadronization for viscous hydrodynamics - in this case the 
Cooper-Frye formula requires a straightforward modification of the form 
of the distribution function /. 

(2) Implementation of more advanced fitting procedures of the femtoscopic 
correlation functions, for instance imaging |149j . We note that even at 
present the output from FEMTO-THERMINATOR contains all necessary in- 
formation for such analyses (the 3-dimensional correlation functions). 

(3) Analysis of the triangular and multipole flow, analysis of the two-particle 
correlations other than those due to the Bose-Einstein statistics. 

(4) Implementation of a more general form (e.g., negative binomial) of the 
multiplicity distributions. This will allow for a realistic modeling of event- 
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by-event multiplicity fluctuations. 

(5) Conservation laws. A modification of the particle-generation algorithm 
might include (approximate) conservation of the global charges and the 
transverse momentum. 

(6) Interface to other formats for the files parameterizing hypersurfaces will 
be developed according to demand of potentially interested users. 
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A Particle data files (particles .data and decays . data) 



The structure of input files containing the information on the particle proper- 
ties and their decays is exactly the same as in the original code [1] — the form 
inherited from SHARE [82J. To maintain integrity, we provide here the most 
basic information. The input file particles . data contains the properties of 
particles such as mass, width, spin, isospin, the quark contents, and the Monte 
Carlo identification number. The file has the format: 

name mass width spin I 13 q s aq as c ac MC 

where: 

name the particle label used in the program, 

mass mass in GeV, 

width width in GeV, 

spin spin, 

I isospin, 

13 3rd component of isospin, 

q, s number of light /strange quarks, 

aq, as number of light /strange antiquarks, 

c, ac number of charmed/anticharmed quarks, 

MC particle's identification number, where available corresponding to the 
standard Monte Carlo particle identification convention [150J. 

The file decays. data contains the information on particle decay channels in 
the format: 

Name parent Name daughtcil Name daughter2 Name daughter3 II] BR C-G?(0/1) 

where BR denotes the branching ratio of the decay and C— G refers to whether 
the branching ratio should be multiplied by a Clebsch-Gordan coefficient 
(0: no, 1: yes). Normally, we use 1 for the two-body decays and for the 
three-body decays, where the full (isospin-dependent) branching ratio BR is 
provided. The entry BR may be used to control the feed-down from the weak 
decays, in particular setting it to would switch off the particular decay chan- 
nel. 



Appears for the three-body decays only. 
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B General definitions 

Variables used to describe the properties of a particle: 



ParticleType.h 


symbol 


deflnition 


mNumber 




particle type number 


mName 




particle name 


mMass 


m 


mass 


mGamma 


r 


width 


mSpin 


s 


spin 


mBarionN 


B 


baryon number 


ml 


I 


isospin 


mI3 


h 


3rd component of isospin 


mStrangN 


S 


strangness 


mCharmN 


c 


charm 


mNq 


n q 


number of u and d quarks 


mNaq 


n q 


number of u and d quarks 


mNs 


n a 


number of s quarks 


mNas 


ri- 


number of s quarks 


mNc 


fle 


number of c quarks 


mNac 


n- 


number of c quarks 


mPDGCode 




particle MC number according to PDG 1 1501 


mMaxIntegrand 




maximum value of the Cooper- Frye integrand 


mMultiplicity 




average multiplity 


mDecayChannelCount2 




number of two-body decay channels 


mDecayChannelCount3 




number of three-body decay channels. 
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Variables used to calculate the Cooper- Frye integrand in the Model_* . cxx 
files. See function GetIntegrand(ParticleType* aPartType). 



in program 


symbol 


definition 


dSigmaP 


p-dY. 


particle's four momentum projected 
on the hypersurface element 


PdotU 


p ■ u 


particle's energy in the fluid local rest frame 


Statistics 


±1 


+1 (for fcrmions), —1 (for bosons) 


space-time 


(Xt,Xx, Xy,Xz) 


(t,x,y,z) 


Cartesian coordinates 


Tau 


T 


proper time 


Rho 


P 


transverse size 


PhiS 




azimuthal angle 


RapS 


Ys 


space-time rapidity 


momentum 


(Pc, Px, Py, Pz) 


(E,p x ,Py,p z ) 


Cartesian coordinates 


Mt 


mr 


transverse mass 


Pt 


Pt 


transverse momentum 


PhiP 


4>p 


transverse-momentum angle 


RapP 


Y 


rapidity 


fluid element 


Ux, Uy 


Ux 1 Uy 


transverse four-velocity 


Vt 


vt 


transverse velocity 


PhiF 


<t>f 


transverse-velocity angle 


RapF 


Yf 


fluid rapidity. 



C Hydro-inspired models 

Common parameters for hydro-inspired models: 



ini file 


in program 


symbol 


description 


RapPRange 


mRapPRange 


yrng 


rapidity range 


RapSRange 


mRapSRange 


v rn s 

1 s 


spatial-rapidity range 


RhoMax 


mRhoMax 


pin ax 


firecylinder transverse size 


Temperature 


mTemp 


T 


freeze-out temperature 


MuB 


mMuB 


PB 


baryon chemical potential 


Mul 


mMul 


PI 3 


isospin chemical potential 


MuS 


mMuS 


PS 


strange chemical potential 


MuC 


mMuC 


PC 


charm chemical potential. 



Integration ranges 
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Y s 



e 

e [o, p m 

e [0, 2tt 



2 ' ~~T~ 



y 6 

Pt £ [0, oo) 
P £ [0, 2?r] 



Ci Cracow Single Freeze- out Model 
Model specific parameters: 



fomodel/krakow.ini 


in program 


symbol 


description 


TauC 


mTauC 


T f 


Cracow proper time. 



Definition of the proper time used in the Cracow model 



Tf = \l t 2 — x 2 — y 2 — z 2 = J t 2 + p 



{C.l) 



C.2 Blast-Wave 

Model specific parameters: 



fomodel/blastwave.ini 


in program 


symbol 


description 


Tail 


mTau 


T 


proper time 


VelT 


mVt 


vt 


transverse velocity of the fluid element. 



C. 3 Blast-Wave A- Class 
Model specific parameters: 



fomodel/bwa.ini 


in program 


symbol 


description 


Tau 


mTau 


r 


proper time 


VelT 


mVt 


vt 


transverse velocity of the fluid clement 


ParA 


mA 


A 


the A parameter in Eq. I|14[l 


Delay 


mDelay 


A 


emission delay 



If the events.ini file has the FreezeOutModel option set to: BWAVLinear, 
BWAVLinearDelay, or BWAVLinearFormation, the value of the fluid transverse 
velocity has a linear form 



vt 



P/Pr, 



(C.2) 



v T + p/p max 

Additionally, for the models BWAVTDelay and BWAVLinearDelay the particle's 
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creation time is shifted randomly according to the exponential distribution 

f(t) = e-^, (C.3) 
where t is the original time of the particle's creation. 



The model BWAVLinearFormation incorporates the particle's formation time 
and shifts the position of the particle, 

Px Pv Pz 

x = x + —t, y = y + -^t, z = z + —t. 



D Freeze-out models based on 2+1 and 3+1 hydrodynamic calcu- 
lations 



Common parameters in hydro-based models: 



lhyquid*.ini 


in program 


symbol 


description 


Rapidity Range 


mRapRange 


yrng 


rapidity range 


*.xml file 

<PARAMETER> 








Tau_i 

Temperature 

Mu.B 

MuJ 

Mu.S 

Mu.C 

device 

colliding_system 

colliding_cnergy 

centrality_min 

centrality_max 

impact_paramctcr 

temperature_at_center 


mTaul 

mTemp 

mMuB 

mMul 

mMuS 

mMuC 

mDeviceName 

mCollidingSystem 

mCollidingEncrgy 

mCentralityMin 

mCentralityMax 

mlmpactParameter 

mTempI 


n 
T 

fJ-B 

M/ 3 
MS 

y/SNN 
Cmax 

b 

Ti 


initial proper time 

freeze-out temperature 

baryon chemical potential 

isospin chemical potential 

strange chemical potential 

charm chemical potential 

accelerator's name, e.g., RHIC, LHC 

colliding nuclei, e.g., Au, Pb 

center-of-mass energy 

centrality lower limit 

ccntrality higher limit 

impact parameter 

central initial temperature 


<VECTOR3D> 








Distance 

DistanceDZeta 

DistanceDPhi 


mDistance 

mDistanccDZcta 

mDistanccDPhi 


d(t,<f>,e) 

dd 
9C 
dd 
8<t> 


distance to the hypcrsurfacc clement 
^-derivative of d 
</>-derivative of d. 



Phase-space integration ranges: 
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2 ' 2 



or e [6» min ,0 max ] 

mm i SmaxJ 
</> G [</> m in, 0max] 



y e 



[0,oo) 
[0,2-tt]. 



Note that in the Lhyquid3D case the space-time rapidity K s is parametrized 
by the 9 angle, F s = (d/A) cos9. 



D.l Lhyqmd 3D 

Model specific parameters: 



*.xml file 

<PARAMETER> 


in program 


symbol 


description 


Lambda 


mLambda 


A 


conversion constant in Y a = (d/A) cos 8 


<VECTOR3D> 








DistanceDTheta 
FluidUx 
Fluidify 
FluidRapidity 


mDistanceDTheta 
mFluidUx 
mFluidUy 
mFluidRapidity 


3d 

ae 

Uy{C,<t>,6) 

Y f (t,<t>,6) 


^-derivative of d 
fluid four-velocity x-component 
fluid four-velocity ^-component 
fluid rapidity. 



D.2 Lhyquid 2D boost-invariant 
Model specific parameters: 



lhyquid2dbi. ini 


in program 


symbol 


description 


RapSRange 


mRapSRange 


y/ ns 


spatial-rapidity range 


*.xml file 

<VECTOR3D> 








FluidVt 
FluidPhi 


mFluidVt 
mFluidPhi 


v T (C,<t>) 


fluid transverse velocity 
fluid velocity azimuthal angle. 
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Table 1 

The contents of the file events . ini 



# 


FreezeOutModel 


i 

1 model parameter file 


# 
# 


KrakowSFO 


_| 

1 . /f omodel/cracow. ini 


# 


BlastWave 


1 . /f omodel/blastwave . ini 


# 


BWAVT 


1 . /f omodel/bwa . ini 


# 


BWAVTDelay 


1 . /f omodel/bwa . ini 


# 


BWAVLinear 


1 . /f omodel/bwa . ini 


# 


BWAVLinearDelay 


1 . /f omodel/bwa. ini 


# 


BWAVLinearFormat ion 


1 . /f omodel/bwa . ini 


# 


Lhyquid3D 


1 . /f omodel/lhyquid3d. ini 


# 


Lhyquid2DBI 


1 . /f omodel/lhyquid2dbi . ini 


# 
# 


Name of the freeze-out model 


# 


available: see table 


above 


# 


default: Lhyquid2DBI 




FreezeOutModel = Lhyquid2DBI 


# 


Custom freeze-out model ini file [not used by default] 


# 


default : 




# 


FreezeOutModellNI = 


/f omodel/lhyquid2dbi . ini 



[Event] 

# Number of events to generate 

# default: 5000 
NumberOfEvents = 5000 



# Event output file format 

# available: root, root&text, text 

# default: root 
EventFileType = root 

[Primordial] 

# Distribution of primordial particles multiplicity 

# available: Poisson 

# default: Poisson 
MultiplicityDistribution = Poisson 

# Number of samples used in determination of primordial multiplicity and max. integrand value 

# default: 5000000 
IntegrateSamples = 5000000 

[Random] 

# Start each event with a new random seed taken from current time (1) or do a constant seed (0) 

# default: 1 
Randomize = 1 

[Directories] 

# Directory with SHARE input files 

# default: ./share/ 
ShareDir = . /share/ 

# Directory with Freeze-Out Model parameter files 

# default: ./fomodels/ 
FreezeOutDir = ./fomodel/ 

# Directory with ROOT macro files *.C 

# default: ./macro/ 
MacroDir = . /macro/ 

# Directory to write the events 

# default: ./events/ 
EventDir = ./events/ 

[Logging] 

# Log file - save information on number of events runed, destination and time 

# default: therminator.log 
LogFile = therminator.log 
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Table 2 

Parameters of f omodel/lhyquid2dbi . ini 
######################################################### 

# Lhyquid 2+1 Boost Invariant # 
######################################################### 

# Rapidity range 

# default: 4.0 
RapPRange =4.0 

# Spatial rapidity range 

# default: 8.0 
RapSRange =8.0 

# location of the hydro-code output containing the definition of the hypersurf ace , 

# velocity profile, thermodynamic properties etc. 

# default: lhyquid2dbi/RHICAuAu200c2030Ti455ti025Tf 145. xml 
FreezeFile = lhyquid2dbi/RHICAuAu200c0005Ti500ti025Tf 145. xml 

# subdirectory to store events of this model. 

# If option is hashed the EventSubDir will be made form FreezeFile name (replace '/' with '-' 

# and remove ".xml") 

# default: lhyquid2dbi-RHICAuAu200c2030Ti455ti025Tf 145/ 
#EventSubDir = lhyquid2dbi-RHICAuAu200c2030Ti455ti025Tf 145/ 



Table 3 

Contents of the femto.ini configuration file. 

[Pair] 

# Type of particle pairs to correlate 

# default: pion-pion, kaon-kaon 
PairType = pion-pion 

[Cuts] 

# Particle time cut [fm] 

# default: 500.0 
TimeCut = 500.0 

[Event] 

# Number of events to mix 

# default: 20 
EventsToMix = 20 

[Switches] 

# Use only primordial particles (yes) or all particles from resonance decays (no) 

# default: no 
EnableOnlyPrimordial = no 

# Enable source histograms 

# default: no 
EnableSourceHistograms = no 

[Logging] 

# Log file - save information on number of events runed, destination and time 

# default: therminator.log 
LogFile = therminator.log 
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Table 4 

Contents of the hbtf it . ini file. 

[Norm_Fit] 

# Normalization initial value 

# default: 1.0 
Norm = 1.0 

# Fit type 

# default: fixed 
NormFitType = fixed 
[Lambda_Fit] 

# Lambda initial value 

# default: 1.0 
Lambda = 1.0 

# Fit type 

# default: free 
LambdaFitType = free 
[Rout_Fit] 

# R_out initial value 

# default : 
Rout =5.0 

# Fit type 

# default: limit 
RoutFitType = limit 

# Rout limits 

# default: 
RoutMin = 2.0 
RoutMax =7.0 
[Rside_Fit] 

# Rside initial value 

# default: 
Rside =5.0 

# Fit type 

# default: limit 
RsideFitType = limit 

# Rout limits 

# default: 
RsideMin =2.0 
RsideMax =7.0 
[Rlong_Fit] 

# Rlong initial value 

# default : 
Rlong =7.0 

# Fit type 

# default: limit 
RlongFitType = limit 

# Rout limits 

# default: 
RlongMin = 4.0 
RlongMax = 14.0 
[QRange] 

# Maximum range of q_out, q_side, q_long 
#default: 0.15 

MaxFitRange =0.15 
[Histograms] 

# Numerator histogram name 

# default: cnuma 
Numerator = cnuma 

# Denominator histogram name 

# default: cdena 
Denominator = cdena 
[Logging] 

# Log file 

# default: therminator.log 
LogFile = therminator.log 



